function [StdParam, StdCovar] = ecmmvnrstd(Data, Design, Covar, Method, CovarFormat)
%ECMMVNRSTD Evaluate standard errors for multivariate normal regression model.
% Given a multivariate normal regression model with missing data, estimate a NUMPARAMS x 1 column
% vector of standard errors StdParam for model parameter estimates and a NUMSERIES x NUMSERIES
% matrix of standard errors StdCovar for covariance parameter estimates, where the model has the
% form
%		Data(k)  ~  N(Design(k) * Param, Covar)
% for samples k = 1, ... , NUMSAMPLES.
%
%	StdParam = ecmmvnrstd(Data, Design, Covar);
%	[StdParam, StdCovar] = ecmmvnrstd(Data, Design, Covar, Method, CovarFormat);
%
% Inputs:
%	Data - NUMSAMPLES x NUMSERIES matrix with NUMSAMPLES samples of a NUMSERIES-dimensional random
%		vector. Missing values are represented as NaNs. Only samples that are entirely NaNs are
%		ignored (to ignore samples with at least one NaN, use MVNRSTD).
%	Design - Either a matrix or a cell-array to handle two distinct model structures. First, if
%		NUMSERIES = 1, Design can be a NUMSAMPLES x NUMPARAMS matrix with known values. This is the
%		"standard" form for regression on a single data series. Alternatively, for any NUMSERIES >=
%		1, Design can be a cell array of either 1 or NUMSAMPLES cells, where each cell contains a
%		NUMSERIES x NUMPARAMS matrix of known values. If Design has a single cell, then it is
%		assumed to be the same Design matrix for each sample. Otherwise, Design must contain
%		individual Design matrices for each sample.
%	Covar - NUMSERIES x NUMSERIES matrix of estimates for the covariance of the residuals of the
%		regression.
%
% Optional Inputs:
%	Method - String to identify method of calculation for information matrix. The choices are:
%		'hessian' - Default method. Use the expected hessian of the observed log-likelihood
%			function. This method is recommended since the resultant standard errors incorporate the
%			increased uncertainties due to missing data.
%		'fisher' - Use the fisher information matrix.
%	CovarFormat - String that specifies the format for the covariance matrix. The choices are:
%		'full' - Default method. Compute the full covariance matrix.
%		'diagonal' - Treat the covariance matrix as a diagonal matrix.
%
% Outputs:
%	StdParam - NUMPARAMS x 1 column vector of standard errors of estimates for each element of the
%		model parameters vector Param.
%	StdCovar - NUMSERIES x NUMSERIES matrix of standard errors of estimates for each element of the
%		covariance matrix Covar.
%
% Notes:
%	WARNING: If calculating standard errors associated with the covariance matrix Covar, i.e., with
%	two output arguments, then this routine is VERY slow.
%
%	Tip - If Method = 'fisher', to obtain more quickly just the standard errors of variance
%	estimates without the standard errors of the covariance estimates, set CovarFormat = 'diagonal'
%	regardless of the form of the covariance matrix.
%
% References:
%	[1] Roderick J. A. Little and Donald B. Rubin, Statistical Analysis with Missing Data, 2nd ed.,
%		John Wiley & Sons, Inc., 2002.
%
%	See also MVNRMLE, ECMMVNRMLE, ECMMVNRSTD.

%	Copyright 2005-2007 The MathWorks, Inc.
%	$Revision: 1.1.6.4 $ $Date: 2007/05/10 13:44:59 $

% Step 1 - check arguments

if nargin < 5 || isempty(CovarFormat)
	CovarFormat = 'full';
else
	if ~any(strcmpi(CovarFormat,{'full','diagonal'}))
		error('Finance:mvnrstd:InvalidCovarianceFormat', ...
			'Invalid format specified for covariance matrix.');
	end
end

if nargin < 4
	Method = 'HESSIAN';
else
	if ~any(strcmpi(Method,{'FISHER','HESSIAN'}))
		Method = 'HESSIAN';
	end
end

if nargin < 3
	error('Finance:ecmmvnrstd:MissingInputArg', ...
		'Missing required input arguments Data, Design, or Covar.');
end

if isempty(Data)
	error('Finance:ecmmvnrstd:EmptyDataArray', ...
		'The required input argument Data is empty.');
end
if isempty(Design)
	error('Finance:ecmmvnrstd:EmptyDesignArray', ...
		'The required input argument Design is empty.');
end
if isempty(Covar)
	error('Finance:ecmmvnrstd:EmptyCovar', ...
		'The required input argument Covar is empty.');
end

[NumSamples, NumSeries, NumParams] = checkmvnrsetup(Data, Design, [], Covar);

if nargout < 2
	MatrixFormat = 'PARAMONLY';
else
	MatrixFormat = 'FULL';
end

% Step 2 - initialization

[CholCovar, CholState] = chol(Covar);
if CholState > 0
	error('Finance:ecmmvnrstd:NonPosDefCovar', ...
		'Covariance matrix is not positive-definite.');
end

% Step 3 - calculate Fisher or Hessian, invert, and pull out std. errors

TestMatrix = ecmmvnrfish(Data, Design, Covar, Method, MatrixFormat, CovarFormat);

Count = sum(~all(isnan(Data),2));

ParamCovar = inv(TestMatrix(1:NumParams,1:NumParams));
ParamCovar = (1.0/Count) .* ParamCovar;

StdParam = sqrt(diag(ParamCovar));

if strcmpi(MatrixFormat,'FULL')
	if strcmpi(CovarFormat,'DIAGONAL')
		CovarCovar = diag(1 ./ diag(TestMatrix(1 + NumParams:end,1 + NumParams:end)));
	else
		CovarCovar = inv(TestMatrix(1 + NumParams:end,1 + NumParams:end));
	end
	CovarCovar = (1.0/Count) .* CovarCovar;

	StdCovar = zeros(NumSeries, NumSeries);

	if strcmpi(CovarFormat,'DIAGONAL')
		for i = 1:NumSeries
			StdCovar(i,i) = sqrt(CovarCovar(i,i));
		end
	else
		k = 0;
		for i = 1:NumSeries
			for j = 1:i
				k = k + 1;
				StdCovar(i,j) = sqrt(CovarCovar(k,k));
				StdCovar(j,i) = StdCovar(i,j);
			end
		end
	end
end
